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We illustrate a technique for fitting lattice QCD correlators to sums of exponentials that is significantly faster 
than traditional fitting methods — 10-40 times faster for the realistic examples we present. Our examples are 
drawn from a recent analysis of the T spectrum, and another recent analysis of the £) — > ;r semileptonic form 
factor. For single correlators, we show how to simplify traditional effective-mass analyses. 
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Most physics results in lattice QCD come from fits of lattice 
correlators to sums of exponentials. For example, we study a 
particular hadron by computing Monte Carlo simulation esti- 
mates G^^{t) of hadronic correlators. 
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with different sources F,, and sinks Tj, that create and de- 
stroy the hadron. The sum over all spatial sites x restricts 
the hadrons to states with zero total three-momentum. Such a 
correlator can be decomposed into contributions from energy 
eigenstates \Ej) in QCD H: 



N 

Gab{t;N) = ^ ajbjexp{-Ejt) 



(2) 



where Ej is the energy, with Ej>Ej^\, and the amplitudes are 
matrix elements, with 



a* = (0|F„(0,0)|£;}, 
fc_,- = (0|F,(0,0)|£;). 



(3) 



The physics is in the energies and the matrix elements, and 
these are determined by fitting fomula dU to the Monte Carlo 
data G^^{t) for a variety of sources and sinks. 

In principle, the number of terms in Eq. is infinite, 
but, in practice, we need only retain a finite number of terms 
because the exponentials suppress high-energy states. The 
number needed depends upon the precision of the simulation 
data G^^, but it is not uncommon to require A^= 10 or more 
terms for good fits to accurate data. The fitting process be- 
comes both cumbersome and time consuming if many corre- 
lators must be fit simultaneously while using such large A^s. In 
this paper we introduce a method that can dramatically sim- 
plify and accelerate such fits. 

The key to this new approach lies in how priors are intro- 
duced. Two types of input data are required for these fits. The 
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first is simulation data for the correlators, consisting of Monte 
Carlo averages G for each a, b and f, and a covariance ma- 
trix that specifies both the statistical uncertainties in each 
average and the correlations between them: 
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{t)^[Gab(t),al,,,y{tf')] 
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This data contributes 

ZlAc{aj.bj,Ej)=Y^ ^ {Gab{t;N)-Gab{t)) 

tM.bt' ,a' ,b' 

0;h,,At,t'){Ga>b'{t';N)-G,,y{t')) 



(5) 



to the function that is minimized by varying parameters a,- 



bj, and Ej in a conventional fit. 



The second type of input data consists of Bayesian priors 
for each fit parameter. Complicated multi-correlator, multi- 
parameter fits are impossible without a priori estimates for 
each fit parameter L2J : 



Ef^Ej±aEy 
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This information is included in a conventional fit by adding 
extra terms to x^{aj,bj,Ej): x^ = Xmc + Xpr where 

Xpvia j,bj,Ej) = 

7=1 I % % ^Ej J 

The priors can also be combined to give a priori estimates for 
the correlators. 



GP^(r;/V)^i:«f^,P'exp(-£fr), 
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where the means and covariance matrix for G^J,(f) are com- 
puted, using standard error propagation, from the means and 
covariance matrix of the priors (Eq. 
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The cost of a traditional analysis goes up rapidly with the 
number of parameters needed to obtain a good fit. In practice, 
however, we are rarely interested in parameters from the large- 
i terms in fit function (|2]i, even when these terms are needed 
for a good fit. Rather than including them in the fit, we can 
incorporate the large- j terms into the Monte Carlo data before 
fitting. To do this, we use the priors to generate an a priori 
estimate for these terms, and then subtract that estimate from 
the Monte Carlo data. This effectively removes the large- j 
terms from the data. Finally we fit the modified data with a 
simpler formula that includes only small-/' terms. 

More explicitly, we can remove terms having n < j<N by 
replacing G^^{t) with (first definition) 



where 



GZ{t;n)^GZ{t)-AGll{t;n), 



j=n+l 



(9) 
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is the j > n part of the fit function. Having removed the 
j >n terms, we can fit G^^{t;n) with the simpler fit func- 
tion, Gah{t',n), rather than Gab{t\N). 

Here we assume that is sufficiently large that AG^^(f;n) 
and therefore &^^{t\n) are independent of to within their 
statistical errors. The covariance matrix for Cj^^{t\n) is 
obtained by adding the covariance matrices of G^^{t) and 
AG^'^(f;n) (that is, adding the errors in quadrature) |4]. 

Removing high- j terms from both the fit function and the 
fit data replaces the original fitting problem — fit an A^-term 
function Gah{t;N) to G^^{t) — by a simpler problem that can 
have far fewer fit parameters: fit an n-term function Gab{t'Ji) 
to G^^(t;n), where n<N. Remarkably, as we showed in yl], 
these two problems are equivalent for high statistics data even 
when n is quite small: that is, fit results (means and stan- 
dard deviations) for the low- j parameters are the same in both 
cases. In the second case, the j > n terms have been "marginal- 
ized," or, in effect, integrated out of the Bayesian probability 
distribution, but in a way that does not affect the analysis of 
the j <n terms. When n<^N, the fit parameters that remain 
are many fewer than what would be required in a standard fit, 
and fitting is much faster 

In this paper we use a variation of this marginalization pro- 
cedure which we find to be more robust when fitting correla- 
tors that fall exponentially quickly with increasing t. In this 
variation the modified correlators are given by (second defini- 
tion) 



Gab (f;n) =G„f, (t) 



G':bit;Ny 



(11) 



which is analogous to the first definition (Eq. (|9]l) but applied 
to the logarithm of the correlator rather than the correlator it- 
self. Again, terms with j>n have been removed, and there- 
fore the modified correlator data can be fit with the simpler fit 
function, Gab{t',n). 
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FIG. 1. Fit X per degree of freedom for sequential fits of 25 T 
correlators with n= 1,2,3 .. . terms in fit function l|2ll. Results are 
plotted versus the cumulative time required for fitting, and are for fits 
of: a) the unmodified simulation data G^^(t) (red circles and dotted 
line); and b) the modified simulation data &^^{t;n) (Eq. (TT)) (blue 
circles and dashed line). The region of good fits is indicated by the 
gray band. 



We now illustrate our new method by applying it to QCD 
simulation data from two recent analyses. For each analy- 
sis, we fit a function, like Gab{t',n), with n terms both to 
untouched simulation data G^^{t), and to modified simula- 
tion data G^^{t;n), from which j > n terms have been re- 
moved using Eq. (fTTT l. We vary «, doing sequential fits with 
n = 1 , 2, 3 . . ., where the best-fit parameter values from one fit 
are used as starting values for the next fit. Sequential fitting 
with increasing « is a standard approach to complicated multi- 
parameter correlator fits; n is increased until the fit's x'^ stops 
changing, at which point enough terms have been include 
to reflect accurately the uncertainties introduced by large- j 
terms. Here we examine the best-fit parameters for each n to 
investigate the rate at which correct results emerge from this 
process. This allows a detailed comparison of our two fitting 
strategies. 

The first data set is a collection of 25 correlators for the 
T(15') meson and its radial excitations (T(25'), T(35), etc.) iU. 
These correlators were made using five different operators for 
both sources and sinks. They were fit to formula (|2]i with pri- 
ors (in lattice units): 

log(£i ) = log(0.3 ± 0. 1 ) = 1 .2 ± 0.3 
log(£;+i =log(0.25± 0.125) = -1.4 ±0.5 



fly = 0.1 ±1.0 



(12) 



except for a local source for which the priors were log(flj) = 
log(0.1 ±0.2) = —2.3 ±2 (local source). These are broad pri- 
ors — more than 100 times broader than the final eiTors for the 
quantities we examine below. We set N ~ 20 when defining 
Gff{t;n) (Eq. (HB); this is roughly twice the size it needs to 
be, but it costs little to make large. In general A^ should be 
chosen so that terms with j>N are negligible compared with 
statistical eiTors. 

In Fig. [T] we plot the per degree of freedom for each 
method versus the time required to get to that value As 
expected, the new algorithm reaches a reasonable x^ with just 
a few terms (n = 2-3), in 20-30 seconds; the traditional al- 
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FIG. 2. Best-fit results from sequential fits of 25 Y correlators with 
n = 1 , 2, 3 . . . terms in fit function (|2}. Results are plotted versus the 
cumulative time required for fitting, and are for fits of: a) the un- 
modified simulation data G^^{t) (red circles and dotted line); and b) 
the modified simulation data G^^{t;n) (Eq. (blue circles and 
dashed line). Results are given for mass splittings between different 
vector 5-states, and for the wave functions at the origin for the lowest 
two states. All results are in lattice units. The gray bands show the 
best-fit result from the modified data after convergence. 



gorithm requires « = 10-11 to obtain a good x^, and 600- 
700 seconds. Similar differences are evident if we look at 
physical quantities extracted from the simulations. In Fig. |2] 
we show results for the 25— 15 mass splitting (in lattice units), 
for the 35—15 mass splitting divided by the 25—15 splitting, 
and for the 15 and 25 mesons' (nonrelativistic) wave func- 
tions at the origin, which come from fit parameters aj for a 
local source. In every case the two algorithms agree on the 
final result, but the new algorithm converges to correct results 
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FIG. 3. Best-fit results from sequential fits of 13 two-point and three- 
point correlators for D and 71 mesons with n= 1,2,3... terms in 
fit function l|2]l. Results are plotted versus the cumulative time re- 
quired for fitting, and are for fits of: a) the unmodified simulation 
data (red circles and dotted line); and b) the modified simulation 
data (Eq. i ll It ) (blue circles and dashed line). Results are given for 
the D-meson mass mo and decay constant fu, and for the k 
scalar form factor at zero recoil momentum /g (0,0,0). All results 
are in lattice units. The gray bands show the best-fit result from the 
modified data after convergence. 



10^0 times faster 

Our second example is from a recent analysis of the D — 
n semileptonic form factor |7]. To extract the form factor 
at four different momenta, this analysis uses a simultaneous 
fit of 13 two-point and three-point correlators: a) a D-meson 
correlator with a pseudoscalar local source and sink; b) four 
TT-meson correlators, one for each pion momentum of inter- 
est, again with local pseudoscalar sources and sinks; and c) 
two three-point correlators D—> /scalar for each of the four 
pion momenta. The fit functions are more complicated for this 
case. For example, the D-meson correlator is fit by a function: 



GD{t;n) 



Y,ajf{Ej,t) 

;=l 



(13) 



where f{Ej,t) = exp(— Eyf) + exp(— ^^(r — t)) is periodic 
with period T — 64, and the second (oscillating) term is due to 
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FIG. 4. The D-meson's effective mass nJeff(?) versus t computed 
from modified simulation data G^*'(/) from whicli every state otlier 
tlian tlie ground state has been removed (using priors). The (very 
thin) gray band shows the weighted average of all /)ieff(?)s, taking ac- 
count of correlations. The thickness of the band indicates the uncer- 
tainty of the average. Note that the largest f s shown here correspond 
to the middle t range. The error bars grow there because meff(f) be- 
comes very sensitive to statistical errors in this region (since periodic 
boundary conditions imply that the derivative of the correlator's non- 
oscillating part vanishes at the midpoint). 



opposite-parity states in the correlator (a feature of staggered- 
quark formalisms like that used in this analysis). The details 
for the other correlators, and the priors are given in [7]. 

Despite the complexity of dealing with both two-point and 
three-point coiTelators, this is a simpler fit than the T case; but 
even here we find that marginalizing most of the fit function 
makes the analysis about 30 times faster We show results in 
Fig.[3]for the D-meson's mass mo and leptonic decay constant 
fo, as well as for the D— > tt scalar form factor /o (0,0,0) at 
zero recoil momentum. All results are in lattice units. Again 
the two approaches agree on the results but the new approach 
has correct results even with only a single term (« = 1) in the 
fit functions. For these fits we set A^= 10 when computing the 
modified data G^f{t;n) (Eq. ( fTTT) ). which is twice as large as 
it needs to be. 

Some insight into how marginalization works can be gained 
by focusing just on the D correlator from this analysis and 
fitting the modified data. 



D V J D \ ) QP^(^t;N) ■ 



(14) 



with only the non-oscillating part of the first term 
in Eq. ( fT3]) — that is, with aif{E\,t). This situation is suf- 
ficiently simple that fitting is not required. The D mass, for 
example, can be obtained by averaging the "effective mass," 

,„..M..rccoshr °"^''t'l:°°^''-'' ). <i5) 



over all f , taking account of correlations between different f s 



The effective mass is plotted as a function of t in Fig. |4] 
It is compared with the weighted average of all 27 meff(f)s 
(gray band), which at rn^^ = 1 . 1584(1 1 ) agrees well with the 
best result, 1.1593(7), from full multi-term fits (top panel in 
Fig.E). 

The first excited state in the D correlator is the opposite- 
parity contribution, which accounts for the oscillation 
in meff(r). Strong statistical correlations between different 
points result in an average meff whose error is more than 
7 times smaller than the best error from an individual meff(f ). 
The eiTors in meff(f) when f < 16 come almost entirely from 
marginalized terms absorbed into the fit data using Eq. ( fT4l l: 
the original Monte Carlo simulation eiTors are negligible 
there. 

Absent marginalization, contributions from excited states 
would limit a traditional effective mass analysis of this data to 
values with t > 16. With marginalization, all f s are used, ex- 
cept for a small number at very small t where the fit function is 
invalid (because of temporal non-locality in the lattice quark 
action). Using 28 fs is possible because we have removed the 
excited states through Eq. ( fT4] i. As a result different meff(f )s 
agree with each other to within their errors: fitting all 27 val- 
ues in Fig.|4]to a constant gives an excellent fit, with a per 
degree of freedom of 0.6. (The result of the fit is, by definition, 
the same as the weighted average reported above.) 

Our new implementation of effective-mass analyses is sim- 
pler and less ambiguous than traditional analyses because we 
are not limited to large fs. More importantly our implemen- 
tation also allows us to quantify the contribution to the uncer- 
tainty in the final ni^^ due to the excited states: here the priors 
for non-oscillating terms in Eq. ( fT3b contribute 0.44<7m, those 
from oscillating terms contribute 0.07(7,„, and the uncertain- 
ties in the Monte Carlo data contribute 0.89(7,,,, where a,„ is 
the standard deviation of m^^^. Such information is essential 
for assessing the reliability of the final result, as well as for 
planning improvements to the analysis. 

In this paper we have shown how to accelerate multi- 
exponential fits to multiple hadronic correlators by removing 
contributions due to excited states from both the fit function 
and the simulation data, before fitting. This technique for 
marginalizing large parts of the fit function greatly reduces the 
number of fit parameters needed in the realistic examples pre- 
sented here, and makes fitting 10-40 times faster Marginal- 
ization also simplifies effective-mass analyses, and general- 
izes easily to analogous multi-state (generalized eigenvalue) 
methods. 
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